germ=read.csv("Establishment.and.survival.csv")

germ$comp = as.factor(germ$comp)
germ$herb = as.factor(germ$herb)
germ$species = as.factor(germ$species)
germ$herb = factor(germ$herb,labels=c("Reduced","Ambient"))
germ$comp = factor(germ$comp,labels=c("Low","Medium","High"))
germ$y = cbind(germ$total.emerged.07.08,400-germ$total.emerged.07.08)
# with(germ,table(comp,herb))

# for compatibility, split into tall and bull datasets
germ.tall = germ[germ$SPECIES=="TT",]
germ.bull = germ[germ$SPECIES=="BT",]

#====================================
#Fit recruitment for Tall Thistle
#====================================
library(lme4)
# fit global model with block effects
germ.tall.mm0 = glmer(y ~ comp + herb + comp:herb+
		(1|block) , data=germ.tall,family=binomial)
germ.tall.mm1 = glmer(y ~ comp + herb + 	comp:herb + 
		(herb|block) , family=binomial, data=germ.tall)
germ.tall.mm2 = glmer(y ~ comp + herb + comp:herb + 	
		(comp|block), family=binomial, data=germ.tall)
	
sapply(list(germ.tall.mm0,germ.tall.mm1,germ.tall.mm2),BIC)
# model 2 is best, so I use (herb|block)

germ.models = list("~comp + herb+herb:comp",	
		"~comp + herb",	
		"~comp",
		"~herb",
		"~1")
germ.tall.fits = lapply(germ.models,function(ff)glmer(paste("y",ff,"+ (herb|block)"),data=germ.tall,family=binomial))
bic = sapply(germ.tall.fits,BIC)		
dbic = bic-min(bic)
wbic = exp(-dbic/2)/sum(exp(-dbic/2))
data.frame(unlist(germ.models),dbic,wbic)[order(dbic),]
germ.tall.fits[[1]]
loglik=sapply(germ.tall.fits,logLik)
k=sapply(germ.tall.fits,function(mm)length(fixef(mm)))+2
stat.res=data.frame(unlist(germ.models),k,loglik,dbic,wbic)[order(dbic),]
write.csv(stat.res,"Germ_res_TT.csv")

# make a plot
nd = data.frame(comp=factor(rep(c("Low","Medium","High"),times=2),levels=levels(germ.tall$comp)),
		herb=factor(rep(c("Reduced","Ambient"),each=3),levels=levels(germ.tall$herb)))

recruit.results = matrix(NA,nrow=nrow(nd),ncol=length(germ.tall.fits))
recruit.se = matrix(NA,nrow=nrow(nd),ncol=length(germ.tall.fits))
for (i in 1:length(germ.tall.fits)){
	# generate predictions from a glmer model using the data from inside the model, mostly. 
	fm = germ.tall.fits[[i]]
	formula = as.formula(germ.models[[i]]) # fixed effects formula
	X = model.matrix(formula,data=nd)
	fe = fixef(germ.tall.fits[[i]])
	vcf = vcov(germ.tall.fits[[i]])

	recruit.results[,i] = X %*% fe
	for (j in 1:nrow(X)){
		pick = which(X[j,]>0)
		recruit.se[j,i] = sqrt(sum(vcf[pick,pick]))
	}

	#results[,i] = 1/(1+exp(-(X %*% fm@fixef)))

}
# model averaged prediction for logbiomass parameter
nd[,"recruitment"] = apply(recruit.results,1,function(x)sum(x*wbic))


t1 = sweep(recruit.results,1,nd[,"recruitment"])
t2 = wbic*(recruit.se^2 + t1^2)
nd[,"recruitment.se"] = apply(t2,1,sum)


# stripplot(total.emerged.07.08~comp|herb,data=germ.tall,jitter.data=TRUE)
#barchart(1/(1+exp(-(recruitment)))~comp|herb,data=nd,ylim=c(0.00,0.25))
# barchart(recruitment~comp|herb,data=nd,ylim=c(0.05,0.14),
# 	panel=function(x,y,...){
# 		panel.barchart(x,1/(1+exp(-y)),...)
# 		pn = which.packet()
# 		pick = nd$herb == levels(nd$herb)[pn[1]] 			
# 		ul = 1/(1+exp(-(y + 1.96*nd[pick,4])))
# 		ll = 1/(1+exp(-(y - 1.96*nd[pick,4])))
# 
# 		panel.segments(x,ll,x,ul,col="black",lwd=2)		
# 	}
# )


#getting the mean and variance of et coeficients, not needed for IPM parameters
coef.names = names(fixef(germ.tall.fits[[1]]))
coef.results = matrix(NA,nrow=length(germ.tall.fits),ncol=length(coef.names))
coef.se = matrix(NA,nrow=length(germ.tall.fits),ncol=length(coef.names))
for (i in 1:length(germ.tall.fits)){
	coef.results[i,] = sapply(coef.names,function(cf)fixef(germ.tall.fits[[i]])[cf],USE.NAMES=FALSE)
	se = sqrt(diag(vcov(germ.tall.fits[[i]])))
	names(se) = names(fixef(germ.tall.fits[[i]]))
	coef.se[i,] = sapply(coef.names,function(cf)se[cf])
}

coef.avg = apply(coef.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(coef.results,2,coef.avg)
t2 = sqrt(coef.se^2 + t1^2)
coef.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

#Parameters for IPM
for (i in 1:6) {IPM.parameters$recruit.ave[IPM.parameters$species=="TT" & IPM.parameters$comp==nd[i,"comp"] & IPM.parameters$herb==nd[i,"herb"]]=nd$recruitment[i]}
for (i in 1:6) {IPM.parameters$recruit.se[IPM.parameters$species=="TT" & IPM.parameters$comp==nd[i,"comp"] & IPM.parameters$herb==nd[i,"herb"]]=nd$recruitment.se[i]}


#====================================
#Fit recruitment for Bull Thistle
#====================================

# fit global model with block effects
germ.bull.mm0 = glmer(y ~ comp + herb + comp:herb+
		(1|block) , data=germ.bull,family=binomial)
germ.bull.mm1 = glmer(y ~ comp + herb + 	comp:herb + 
		(herb|block) , family=binomial, data=germ.bull)
germ.bull.mm2 = glmer(y ~ comp + herb + comp:herb + 	
		(comp|block), family=binomial, data=germ.bull)
	
sapply(list(germ.bull.mm0,germ.bull.mm1,germ.bull.mm2),BIC)
# mm1 is best, so I use (herb|block)

germ.models = list("~comp + herb+herb:comp",	
		"~comp + herb",	
		"~comp",
		"~herb",
		"~1")
germ.bull.fits = lapply(germ.models,function(ff)glmer(paste("y",ff,"+ (herb|block)"),data=germ.bull,family=binomial))
bic = sapply(germ.bull.fits,BIC		)
dbic = bic-min(bic)
wbic = exp(-dbic/2)/sum(exp(-dbic/2))
data.frame(unlist(germ.models),dbic,wbic)[order(dbic),]
germ.bull.fits[[3]]
loglik=sapply(germ.bull.fits,logLik)
k=sapply(germ.bull.fits,function(mm)length(fixef(mm)))+2
stat.res=data.frame(unlist(germ.models),k,loglik,dbic,wbic)[order(dbic),]
write.csv(stat.res,"germ_res_BT.csv")


# make a plot
nd = data.frame(comp=factor(rep(c("Low","Medium","High"),times=2),levels=levels(germ.bull$comp)),
		herb=factor(rep(c("Reduced","Ambient"),each=3),levels=levels(germ.bull$herb)))

recruit.results = matrix(NA,nrow=nrow(nd),ncol=length(germ.bull.fits))
recruit.se = matrix(NA,nrow=nrow(nd),ncol=length(germ.bull.fits))
for (i in 1:length(germ.bull.fits)){
	# generate predictions from a glmer model using the data from inside the model, mostly. 
	fm = germ.bull.fits[[i]]
	formula = as.formula(germ.models[[i]]) # fixed effects formula
	X = model.matrix(formula,data=nd)
	fe = fixef(germ.bull.fits[[i]])
	vcf = vcov(germ.bull.fits[[i]])

	recruit.results[,i] = X %*% fe
	for (j in 1:nrow(X)){
		pick = which(X[j,]>0)
		recruit.se[j,i] = sqrt(sum(vcf[pick,pick]))
	}
	#results[,i] = 1/(1+exp(-(X %*% fm@fixef)))
}
# model averaged prediction for logbiomass parameter
nd[,"recruitment"] = apply(recruit.results,1,function(x)sum(x*wbic))

t1 = sweep(recruit.results,1,nd[,"recruitment"])
t2 = wbic*(recruit.se^2 + t1^2)
nd[,"recruitment.se"] = apply(t2,1,sum)

# stripplot(total.emerged.07.08~comp|herb,data=germ.bull,jitter.data=TRUE)
# barchart(1/(1+exp(-(recruitment)))~comp|herb,data=nd,ylim=c(0.00,0.25))


#getting the mean and variance of et coeficients, not needed for IPM parameters
coef.names = names(fixef(germ.bull.fits[[1]]))
coef.results = matrix(NA,nrow=length(germ.bull.fits),ncol=length(coef.names))
coef.se = matrix(NA,nrow=length(germ.bull.fits),ncol=length(coef.names))
for (i in 1:length(germ.bull.fits)){
	coef.results[i,] = sapply(coef.names,function(cf)fixef(germ.bull.fits[[i]])[cf],USE.NAMES=FALSE)
	se = sqrt(diag(vcov(germ.bull.fits[[i]])))
	names(se) = names(fixef(germ.bull.fits[[i]]))
	coef.se[i,] = sapply(coef.names,function(cf)se[cf])
}

coef.avg = apply(coef.results,2,function(x)sum(x*wbic,na.rm=TRUE))
t1 = sweep(coef.results,2,coef.avg)
t2 = sqrt(coef.se^2 + t1^2)
coef.avg.se = apply(t2,2,function(x)sum(wbic*x,na.rm=TRUE))

#Parameters for IPM
for (i in 1:6) {IPM.parameters$recruit.ave[IPM.parameters$species=="BT" & IPM.parameters$comp==nd[i,"comp"] & IPM.parameters$herb==nd[i,"herb"]]=nd$recruitment[i]}
for (i in 1:6) {IPM.parameters$recruit.se[IPM.parameters$species=="BT" & IPM.parameters$comp==nd[i,"comp"] & IPM.parameters$herb==nd[i,"herb"]]=nd$recruitment.se[i]}